TRAJECTORY OPTIMIZATION 
USING EAST -TIME REPETITIVE COMPUTATION 
James S. Raby 

Ames Research. Center, NASA, Moffett Field, Calif. 

and 

Rodney C. Wingrove 

Ames Research Center, NASA, Moffett Field, Calif. 

ABSTRACT 


A perturbation method of computing optimum trajectories is described. 
This method uses fast -time repetitive computations in determining control 
impulse response functions and requires only dynamic solutions of the 
state equations. The solution of additional linear adjoint equations 
is not required. 

Both a hybrid and a digital computer mechanization of this impulse 
response method are described. Sample solutions for each computer mecha- 
nization are presented, using the same problem formulation in each case. 
The problem constraints -were that a reentry vehicle travel a specified 
range, that the control remain within specified limits, and that the 
heat input to the vehicle be minimized. A comparison is made between 
the hybrid computer and the digital computer mechanizations, each computed 
near optimum trajectories in about 2 minutes and full optimum trajec- 


tories in about 5 minutes of computer time. 
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INTRODUCTION 

Space vehicle trajectories must be near optimum in the sense that 
some parameter is either a maximum or a minimum; for example, in reentry 
the trajectory to desired terminal conditions is near optimum when the 
total aerodynamic heating is a minimum. Several perturbation methods, 1 
such as the calculus of variations, applications of the maximum principle, 
and direct steepest descent, have been considered for determining the 
time histories of nonlinear controls that correspond to optimum trajec- 
tories. 

The computations 2 in these previous optimization studies involved 
the dynamic solution of two sets of equations : (l) nonlinear state 

equations and (2) linear adjoint equations. An alternate perturbation 
computation technique - the impulse response method 3 - will be discussed 
here. This method differs from previous studies in that only the solu- 
tion of the nonlinear state equations is used. The response of given 
functions (e.g., terminal error or quantity to be optimized) to a con- 
trol impulse is determined along the trajectory by fast -time repetitive 
computations rather than by a solution of adjoint equations. This 
impulse response method enables the investigator to retain an intuitive 
understanding of the optimization process. Furthermore, since adjoint 
equations are not required, the state equations or cost functions need 
not be amenable to linearization. The impulse response method does 
require many solutions of the state equations; however, the programming 
is straightforward and the task of computing a large number of dynamic 
solutions is ideally suited to modern high-speed computers. 
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NOTATION 

The following notation is used in the body of the text. Additional 
symbols are described as they are introduced. 

L/D control value of lift -drag ratio 
n number of storage points in control time history 

t time 

tf final time 

t 0 initial time 

At time increment of control impulse 

u control 

Au height of control impulse 

cp cost at final time 

Ap(t) change in cost at final time due to control impulses at time t 
\|r state value at final time 

desired state value at final time 

At(t) change in state value at final time due to control impulses at 
time t 

GENERAL OUTLINE OF METHOD 

The impulse response method, as discussed in this report, uses the 
steepest descent optimization process . 4 " 7 The process commences with 
any nonoptimal trajectory from which a slightly improved one is derived. 
The improved trajectory is then used as a new nominal trajectory, and 
the procedure is repeated until the optimum or nearly optimum trajectory 


is found. 
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The iterative procedure is: (l) estimate a reasonable nominal 

control program; (2) determine impulse response functions that indicate 
the best method of making small changes in the control that will decrease 
the cost (the quantity to be minimized); (3) compute a new nominal con- 
trol by adding this change in control to the previous nominal control 
(this results in a new trajectory with a decreased cost) ; ( 4 ) repeat 
step 2. This iterative process continues until the change in cost for 
each new trajectory is very small; the control is then very near a 
local optimum. If at any point along the trajectory a limit value of 
the control is reached before the cost is completely minimized, no 
further optimization is possible at that point. In this case, the pro- 
cess continues until at each point on the trajectory either a local 
optimum or the control limit is reached. 

Computation of Impulse Response Functions 
The technique by which the impulse response function is determined 
is the most important feature of the impulse response method. Figure 1 
illustrates the manner in which the influence of small control changes 
on the cost are calculated. The equations of motion are first solved 
with a positive control impulse at time t superimposed upon the nominal 
control. During the next solution of the equations of motion, a nega- 
tive control impulse of the same magnitude is inserted at time t . The 
impulse response, Ap, is derived from these two solutions. In a similar 
manner, the impulse response is determined at successive times along 
the trajectory. The impulse response function, Ap(t) , is the complete 
time history of Ap. Its computation for the same control impulse at 
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different times along the trajectory is defined as one iteration. 

This corresponds to previous optimization studies 4-7 that vised one itera- 
tion of the adjoint equations to compute essentially the same impulse 
response function along the trajectory. 


Calculation of Minimum Cost 

When the cost is to he minimized and there is no terminal constraint, 
the impulse response function is used in the steepest descent technique 
to modify the control toward the optimum in the following manner : 

( New \ / Previous\ 


Nominal 
\ Control/ 


Nominal 
\ Control / 


K<p Ap(t) 


( 1 ) 


The gain K^> weights the impulse response function for the cost; 
its sign is negative to decrease the cost. The magnitude of is deter- 
mined experimentally for each problem: too large a gain may cause 

instability in the convergence procedure, while too small a gain may 
extend the time of convergence. 

A representative sample of what one may expect with this type of 
optimization procedure is sketched in figure 2. During the first itera- < s Pig. 2 
tion, the repetitive solutions determine the impulse response function, 

Acp(t) . This Acp(t) is added to the nominal control with an appropriate 
gain K^, and the new nominal control time history, as shown in the center 
of figure 2, is obtained. The iterative process is repeated until the 
optimum control is reached. The optimum control may take on either or 
both of the properties illustrated in the final iteration of figure 2. 

In the region (A) the impulse response function is, for all practical 
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purposes, zero. This implies that a small change in control in this 
region will not modify the cost; thus the control is at a local optimum. 
In region (B) the control is at the limiting constraints, and the 
impulse response function indicates that only control beyond the con- 
straint will decrease the cost. Thus, on the constraint, the control 
is at a local optimum. 


Minimum Cost With Terminal Constraint 
When a terminal constraint (or destination) is to be reached, 
while minimizing the cost, the iteration procedure is performed as 
follows : 


{ New 
Nominal 
Controls 


V 


Previous^ 
Nominal 
V Control / 


+ K^p Acp(t) + Kji, Av(t) 


( 2 ) 


where denotes the state variable at the final time. The quantity 
Aty(t) represents the change in the state variable due to control impulses, 


and is evaluated in the same manner as Ao(t) . 

Gains Kq> and are constants for each Iteration. Gain K<p 
weights the impulse response function for cost; its sign is negative 
to decrease the cost. Gain must be calculated for each iteration 

so that the term Aty(t) will account for terminal displacement due to 
the optimizing term, Ky Acp(t) , and correct any terminal displacement 
error from the previous iteration. 
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The equation for is: 

Pf 

j* Ap(t)At(t)dt 

to % - +• 

K= -^p— + 2 Au At — (3) 

J f At 2 (t)at J f Aij/ 2 (t)dt 



Steepest descent Terminal error 

optimization term correction term 


The derivation of this equation can be found in appendix A. Au 
is the height of each control impulse; At is the time interval of 
each control impulse; ty* is the desired end-point value; and i|/q - \]/ 
represents any terminal displacement error from each previous itera- 
tion. This equation gives the general form of the steepest descent 
computations; the computer mechanization of this method will be dis- 
cussed next. 

COMPUTER MECHANIZATION AND RESULTS 

The impulse response method has been mechanized on both a hybrid 
and a digital computer to determine the optimal time history of the 
lift -drag ratio (control L/D) that must be flown for a vehicle return- 
ing into the earth’s atmosphere. The example problem requires that 

(1) the cost, cp, which is the heat input to the vehicle, be minimized, 

(2) the vehicle arrive at a terminal constraint, (destination) , 
and (3) the control time history remain within specified limits. The 
solution to this particular problem is known a priori to be a bang- 
bang control; therefore, the final results can be verified. 

The equations of motion are presented in appendix B. The vehicle 
characteristics and flight conditions were those of a manned capsule 
returning from earth orbit and having the following parameters: 
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Initial conditions: altitude 250,000 ft 

horizontal velocity 25,000 fps 

vertical velocity 748 fps 

range to destination 1,000 miles 

Stopping conditions: altitude 100,000 ft 

Control limits: L/D 0 < L/D < 0.5 

The computer systems used in the two mechanizations are described 
in appendix C. 

Hybrid Computer Mechanization 

The major elements of the hybrid computer consisted of: (l) an 

analog computer to solve the trajectory equations, (2) parallel digi- 
tal logic units to control the computer program, (3) delay line memories 
to store the control time history, and (4) D-A and A-D converters to 
transfer the control time history between the analog computer and the 
delay line memory. 

The L/D time history was stored in 64 word serial delay line 
memories with a resolution of 13 bits. The access time of the serial 
memory was 128 psec. To permit a complete solution of the trajectory 
equations within the 128 psec, the analog computer was time -scaled at 
3750 to 1. 

The mechanization of this problem on the hybrid computer is 
illustrated in figures 3 and 4: Figure 3 is the problem flow chart, 

and figure 4 illustrates the logic used in controlling the problem. 

The serial memory unit is continuously driven by counter pulses (Logic 
No. 1), The output of the serial memory is the nominal control time 
history with n points. This time history is used, together with the 
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appropriate control impulse, to solve the trajectory equations. These 
equations are started at the specified initial conditions with Logic 
No. 2, and stopped with Logic No. 3 when the trajectory reaches the 
specified end condition on altitude. The final values of the cost 
quantity (heat) and the state quantity (range) are stored at the 
end of each run as indicated by Logic Nos. 4 and 5* The positive or 
negative control impulse is added to the nominal control input with 
Logie Nos. 6 and J , respectively. Logic No. 8 inserts the modifying 
control (K<p Ap(t) + At(t) ) into the serial memory. This procedure 
runs in essentially a continuous manner, that is, one point out of the 
n points in the nominal control history is updated after each two repet- 
itive computations. After 2n repetitive computations (one iteration), 
every point in storage has been modified and the process is repeated. 

For each iteration, gains K<p and are held constant. As previously 
mentioned, gain determines the relative speed and stability of'- the 

convergence onto the optimum. The corresponding value of to be used 

with each new iteration is calculated by equation (3) as a function of 
the terminal error from each previous iteration (t^ - t) and the follow- 
ing two integrated quantities from each previous iteration: 

tf 

J’ Acp(t)At(t)dt (4) 

to 

and 

tf 

J At (t)dt 


( 5 ) 
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Time t Q was represented "by a logic signal at the first repetitive 
computation in an iteration cycle and time tf was represented hy a 
logic signal at the last computation in an iteration cycle. It should 
be noted that during those parts of the trajectory when the control was 
at a constraint limit, no further optimization was possible and the 
integration of equations (4) and (5) was not carried out during those 
times . 

Hybrid Computer Results 

The results obtained from the hybrid simulation are illustrated in 
figures 5 and 6. Figure 5 shows a portion of one iteration, while fig- 
ure 6 shows the convergence to the optimum control L/D. 

In the upper trace of figure 5? the control impulses are super- 
imposed upon the initial nominal control. Each control impulse had a 
magnitude of L/D = ±0.25 and a time increment of one clock pulse 
(0.002 sec). This control impulse was chosen because it gave variation 
in the final range and heat load on the order of ±5 percent. The inte- 
grated heat loads along each of the repetitive trajectories are presented 
in the next trace. The difference between the final quantities for each 
pair of subsequent runs is Ap, and represents the heat load impulse 
response. 

In figure 6, the first few iterations of the converging optimization 
procedure are illustrated together with the thirtieth iteration. In the 
upper trace the nominal control is recorded as it is read out of serial 
memory every 128 +1 counter pulse (with Logic No. 8). This gives a con- 
venient time history to show the manner in which the control has been mod 
ified during each iteration. Notice that the control is limited within 
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0 < L/D <0.5- This was achieved by simply limiting the output of the 
serial memory to these values. The modifying control shown in the lower 
trace of figure 6 is the sum, Ap(t) + Ify Aty(t) . For this series of 
runs, a constant = -2 . 5x10 ~ 3 / Bt u/ ft 2 permitted fairly rapid conver- 
gence while program stability was maintained. The value of was 

calculated for each iteration by equation (3) to be that value which 
kept the final value of range near 1,000 miles. 

As can be seen in figure 6, the optimum control variation for this 
particular example was a bang -bang control. With the steepest descent 
method, it was found that near -optimum control could be achieved in the 
first few iterations, but to "square up the corner" and achieve full 
optimum control required more iterations (20 to 30.) • 

Digital Computer Mechanization 

The major elements of the digital computer system consisted of: 

(l) a digital computer to solve the trajectory equations, perform the 
logical control of the program, and store the control L/D time histo- 
ries; (2) a line printer to print hard copies of the results; and 
(3) D-A converters and a strip chart recorder for fast observation of 
trends . 

The digital program was written in floating point symbolic language. 
Since the optimization technique requires repetitive computation of the 
trajectory, the choice of an integration routine was very important. A 
fast, stable, and fairly accurate routine was needed. These requirements 
conflict to some extent; 8,9 however, the fourth-order Adams -Bashford inte- 
gration algorithm gave satisfactory results at a step size of 5 seconds, 



J. S. Raby -Trajectory Optimization 


12 


provided a satisfactory starter was used. The starter used the lower 
order Adams -Bashf or d algorithms with a step size of 1 second. 

The program flow is as follows (see fig. 7); (l) A nominal control \Pig. 7 

time history is used to calculate the nominal trajectory. (2) This tra- 
jectory is stored for use as the initial conditions for the repetitive 
computations of the trajectory. (3) At the initial point along the nom- 
inal trajectory, the control is perturbed with a positive pulse, and a 
new trajectory is calculated. (!)•) At this same point on the trajectory, 
the control is perturbed with a negative pulse and another new trajectory 
calculated. (5) From these two repetitive computations of the trajectory 
the heat impulse response, Acp, and the range impulse response, Aty, are 
calculated. (6) The program is then advanced to new initial conditions 
along the nominal trajectory by the length of the integration step size. 

(7) Steps (3) through (6) are repeated until the initial altitude reaches 
the stopping condition (100,000 ft). (8) At this time, a new nominal 
control time history is computed using equations (2) and (3) • (9) Steps 

(l) through (8) are repeated. This iterative computation continues until 
an optimum trajectory is reached. 

Digital Computer Results 

The first five iterations and the twentieth iteration of the 
digital simulation are illustrated in figure 8. The upper trace of < \Fig . 8 

figure 8 shows the control L/D time history. During the first itera- 
tion, the control L/D was a constant 0.25; at the end of this itera- 
tion it was modified by equations (2) and (3) • By the fifth iteration 
the control L/D was approaching bang -bang and by the twentieth iteration it 


was essentially bang -bang. The pulse used to perturb the trajectory 
had a height of 0.25 L/D and a width equal to one integration 
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step size. For this pulse, a constant value of IQp = -7 • 5X10 -2 / Btu/ ft 2 
permitted a fairly rapid convergence and the computation remained quite 
stable . 

The second trace of figure 8 shows the variation in heat from one 
iteration to the next. The heat, which is the cost in this example, 
decreases markedly during the first five iterations and nearly reaches its 
final value by the end of the fifth iteration. Tables I, II, and III give 
the results in tabular form. The range is shown to remain near 1,000 miles 
while the heat is reduced from 23,491 Btu/ft 2 , at the end of iteration 1, 
to 21,517 Btu/ft 2 at the end of iteration 5* The major change during iter- 
ations 5 through 20 was to "square up" the L/D control and achieve the 
full optimum control. At the end of iteration 20, the final range achieved 
was 999*9 miles and the heat 20,966 Btu/ft 2 . During the optimization pro- 
cedure the range varied slightly about the desired value of 1,000 miles 
and the heat load was reduced about 10 percent. 

Discussion of Hybrid and Digital Results 

It was interesting to observe that both the hybrid and the digital 
simulations required approximately the same amount of computer time, 
approximately 2 minutes to obtain near optimum trajectories and approx- 
imately 5 minutes to obtain full optimum trajectories. However, it 
should be pointed out that no real attempt was made to minimize either 
of these computing times . There are several methods for reducing the 
computer time required to obtain optimum trajectories. One method would 
be to select the gain K<p automatically for each iteration instead of 
using a constant value for the entire computing run. This would cause 
the solution to converge to an optimum in fewer iterations at the expense 
of complicating the computer program. Another method of decreasing 
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the computation time would be to decrease the number of points used to 
store the control time history which would decrease the number of repet- 
itive computations required for each iteration. 

The results obtained by both the hybrid and the digital computer 
appear satisfactory for engineering purposes. The final values of range 
and heat computed by the two simulations agree to within approximately 
1 percent and both simulations arrived at the same bang -bang control time 
histories . 

One excellent feature of the digital simulation was the program 
documentation obtained by using the on-line typewriter and line printer. 

The typewriter documented every change made during the time the program 
was in the computer, and the line printer permitted the analysis of each 
variable at specific points along the trajectory. Equally valuable was 
the strip chart recording normally obtained in hybrid computation. It 
was obtained in the digital program by D-A conversion of the digital 
variables. This "quick look" capability made it possible to observe 
trends not readily apparent in numerical printouts. 

The result of this test example was no surprise. In simulations 
that require complicated logic control of the program and a moderate amount 
of storage, there is a distinct advantage. to using a digital computer. 

It proved reasonable to use a digital computer in this simulation because 
there was only a moderate number of sing>lified equations to be solved. 

If the number of equations were increased, the time to solve them on the 
digital computer would, of course, also increase. 
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COMPARISON WITH ADJOINT STEEPEST DESCENT 

A current reentry optimization study at Ames Research Center is using 
both the impulse response method of this report and the standard adjoint 
steepest descent computing method. This study is of interest because the 
two methods have been programmed on the same computer (IBM 709^) and their 
ability to solve several identical problems has been compared. 

Representative solutions obtained from the two methods are illustrated 

in figure 9* This particular example is for the same reentry vehicle and 

initial flight conditions used in the previous example of this report. 

However , the cost function is of the form: 

tf 

cp = J [(Heat rate) + (Hrag) 2 ]dt 
to 

and there is no terminal constraint. This was chosen in order to illus- 
trate a problem formulation that does not represent a bang -bang optimal 
control result. 

The results of the twentieth iteration are shown in figure 9* The 
upper curve shows that the control solutions are almost identical. In 
the lower curve the impulse response function Acp(t) has been normalized 3 
for comparison with the corresponding results obtained by the adjoint 
solutions. Figure 9 demonstrates that the two methods arrive at essentially 
the same final solution. 

For reentry problems similar to the one presented herein, it has been 
found that the computing time required with the adjoint method is about 
one order of magnitude less than that required by the impulse response 
method. Because the adjoint method uses less computer time, it has been 
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the more desirable method for production runs that require a large number 
of optimized trajectories. However , because the impulse method is straight- 
forward to program and because the engineer is able to retain an intuitive 
understanding of the optimization procedure, the impulse method has been 
the more desirable method for initial problem mechanization. Furthermore, 
adjoint equations require linearization and, therefore, cannot be used in 
some problem formulations. For example, in reentry problem formulations 
with complicated heat -balance equations , 10 rather than the simple heating 
expression shown in appendix B, the heat rate cannot be linearized. In 
this type of formulation, the impulse response method has provided the 
only practical solution.* 

*lynamic programming was also tried for this problem but the computer 
time was found to be excessive, one to two orders of magnitude greater 
than that required with the impulse response method. 

CONCLUSIONS 

This paper has described reentry trajectory optimization using the 
impulse response method. The method requires that the computer perform 
a large number of fast -time repetitive computations in solving the state 
equations and in determining impulse response functions. These repeti- 
tive computations are readily performed by both hybrid and digital 
computers . 

The mechanization of the inpulse response method on both hybrid and 
digital computers was found to be straightforward. Near optimum reentry 
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trajectories were obtained in approximately 2 minutes and full optimum 
reentry trajectories in approximately 5 minutes of computer time. The 
solutions obtained from either mechanization agreed to within approximately 
1 percent. 

The impulse response method has been compared with the adjoint steepest 
descent method. The solutions obtained by either method were essentially 
identical. The adjoint method requires less computer time; however, the 
impulse response method does not require familiarization with or -use of 
an auxiliary set of linear adjoint equations. Furthermore, for problem 
formulations that axe not amenable to linearization, the impulse response 
method may be the only practical method. 
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APPENDIX A 

DERIVATION OP EQUATION FOR !% 

Along a normal trajectory , small changes, St, in the terminal state 
due to small changes, Su(t) , in control can be approximated by: 


St 


1 

2 Au At 



Su(t)At(t)dt 


(Al) 


where Au is the height of each control impulse and At is the time 
interval of each control impulse. Substituting IQp Ap(t) + At(t) 
from equation (l) for Su(t) , we have: 


St = b J tf [K 9 + % A^Ct) ]dt (A2) 

to 

Solving for and letting -St - t d - t (the previous terminal error) , 

we obtain: tf 

J Acp(t)At(t)dt 

t 0 t d - t 

— +2AuAt lI^ (A3) 

J At 2 (t)dt J At 2 (t)dt 

v t 0 t 0 

v J v — — ■ 

Steepest descent Terminal error 

optimization term correction term 
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APPENDIX B 

REENTRY TRAJECTORY EQUATIONS 


The following simplified equations derived for flight within the 
atmosphere were used for the example problem herein. The primary assump- 
tions include a spherical nonrotating earth, small flight -path angles, and 
a constant gravity term. The derivation of these equations and their 
applicability have been considered in a number of reports. 11 


•• V 2 

h = -g + — + 


'CpA' 

m 


I 


where 


t = -(tt) I ^ 


/ 2 


[.tf 

\|/ = J V dt 
to 

tf 

cp = 1.7x10 “ 8 J s/^ V 3 dt 

to 


CpA. 

m 

g 

h 

L 

D 

r 

V 

P 

9 

t 


drag loading, 2.0 ft 2 /slug 

local gravitational acceleration, 32.2 ft/sec 2 
altitude, ft 

control value of lift -drag ratio 

radius from earth center, 21.1x10 s ft 
horizontal velocity, fps 

atmosphere density, 0,00237 e “k/ 23 > 500 slug/ft 2 
total heat input, Btu/ft 2 
final range, ft 
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APPENDIX C 

DESCRIPTION OP COMPUTER SYSTEMS 

In order to make meaningful a comparison of the results obtained 
from the analog and digital simulations, it is necessary to very briefly 
describe the computer systems used. 

The analog computer was an EAI 231R-V equipped with electronic 
mode control of the amplifiers. The logic element of the hybrid simu- 
lation was an EAI DOS 350. The DOS 350 has a patchboard which permits 
one to combine logical elements, such as AND gates, flip-flops, shift 
registers, counters, etc., into complicated logic systems. It also has 
several delay line memories of various lengths as well as A-D and D-A 
converters for communicating between the DOS 350 and the analog computers. 

The digital computer was an EAI 8400 mode 0 computer which had a 
2 psec memory access time, an average floating point add time of approx- 
imately 13 psec, an average floating point multiply time of approximately 
15 psec, and a floating point word size of 32 bits. 
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TABLE I. - ALTITUDE TIME HISTORIES TABLE II. - CONTROL TIME HISTORIES 


Time, 

sec 


Altitudes 

10 3 ft 


Time, 

sec 


Control 

l/d 


1 

iteration 
2 10 

20 

1 

Iteration 
2 10 

20 

0 

250 

250 

250 

250 

"o' 

0,250 

0.212. 

0 

b 

6o 

207 

206 

197 

197 

60 

.250 

.192 

0 

0 

120 

180 

184 

160 

160 

120 

.250 

.291 

.500 

.500 

180 

182 

184- 

212 

209 

180 

.250 

.290 

.500 

.500 

240 

158 

165 

196 

195 

240 

.250 

.294 

.447 

.500 

300 

125 

132 

142 

142 

300 

.250 

.294 

.500 

.500 

360 



127 

127 

36O 


.250 

•347 

.475 

400 



104 

110 

400 



.254 

.277 


TABLE III. - TERMINAL CONDITIONS 


Iteration 



' 1 

2 

5 

10 

20 

Time, sec 

344 

358 

389 

407 

4l4 

Altitude, 

£ 

CO 

99.3 

99-9 

99.8 

99.5 

98.8 

Range, miles 

997-7 

1001. 5 

1003.7 

1002.8 

999.9 

Heat, 

Btu/ft 2 

23491 

23197 

21517 

21025 

20966 
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FIGURE CAPTIONS 

Figure 1.- Computation of the impulse response. 

Figure 2.- Computation of the optimal control. 

Figure 3 • - Hybrid computer flow diagram. 

Figure 4. - Hybrid computer program logic. 

Figure 5 • - Hybrid repetitive computations . 

Figure 6.- Hybrid computation of the optimal control. 

Figure 7 .- Digital computer flow chart. 

Figure 8.- Digital computation of the optimal- control. 

Figure 9*- Comparison of the impulse response and adjoint steepest 
descent methods. 
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